#install.packages("xlsx", "plyr", "data.table", "ggplot2", "RColorBrewer", "stringr", "ggpubr")

library(xlsx)
library(plyr)
library(data.table)
library(ggplot2)
library(RColorBrewer)
library(stringr)
library(ggpubr)

setwd("E://OneDrive//研究//8-学位论文//数据分析处理//第1章")  #设置数据读取与输出文件夹

##################
###  数据输入  ###
##################
rm(list=ls())
Data <- read.xlsx("1-data.xlsx",sheetName = "data", encoding = "UTF-8")  #读取待分析数据
M_ECFs <- data.frame(m = c("DC", "DN", "NC", "TC", "TD", "TN", "PN"),
                     s1 = c("D", "D", "N", "T", "T", "T", "P"),
                     s2 = c("C", "N", "C", "C", "D", "N", "N"))  #两ECFs组合关系

############################
###  创建自定义统计函数  ###
############################

bootstrap <- function(data, tv, cv, group_c = "Target_Type", group_t, n_iter, n_samp) {  
  #data：需数据框格式, tv：处理变量, cv：对照变量, group_c：对照分组（默认物种）, group_t：处理分组, n_iter：迭代次数, n_samp：抽样数
  library(plyr)                                #dlply拆分数据框
  library(data.table)                          #data.table包tstrsplit、rbindlist函数
  
  group <- c(group_t, group_c)                 #分组变量
  
  data_c <- dlply(data, group_c)               #按组将数据拆分进list
  C_names <- names(data_c)                     #将分组组合
  C_g <- as.data.frame(tstrsplit(C_names, ".", fixed=TRUE))
  #分组组合拆为数据框
  names(C_g) <- group_c                        #命名分组数据框列名
  
  
  data_t <- dlply(data, group)                 #按组将数据拆分进list
  T_names <- names(data_t)                     #将分组组合
  T_g <- as.data.frame(tstrsplit(T_names, ".", fixed=TRUE))
  #分组组合拆为数据框
  names(T_g) <- group                          #命名分组数据框列名
  
  
  E <- list()                                  #效应list
  for (i in c(1:length(data_t))) {                 #按分组数循环
    bs <- c()
    for (j in c(1:n_iter)) {
      set.seed(j)
      bs <- append(bs, mean(sample(data_t[[i]][tv][,1], n_samp, replace = TRUE), na.rm = T) - mean(sample(data_t[[i]][cv][,1], n_samp, replace = TRUE), na.rm = T))
    }
    E[[i]] <- cbind(data.frame(T_g[i,]), bs, row.names = NULL)
    names(E[[i]]) <- c(group, "E")             #给数据框命名
  }
  return(E)
}

mq <- function(data, v, group) {               #将data数据框按group分组计算v的平均值mean，2.5%和97.5%的分位数down、up，和大于0的比例p
  library(plyr)                                #dlply拆分数据框
  datac <- ddply(data, group,
                 .fun = function(xx, col) {
                   c(mean = mean     (xx[[col]], na.rm=T),
                     down = quantile (xx[[col]], 0.025, na.rm=T),
                     up   = quantile (xx[[col]], 0.975, na.rm=T),
                     p    = length(which(xx[[col]]>0))/length(xx[[col]])
                   )
                 },
                 v
  )
  names(datac) <- c(group, "mean", "down", "up", "p")
  datac <- datac[which(is.nan(datac$mean) == F),]
  return(datac)
}

replace <- function(vector, from, to) {       #批量替换向量vector中的元素，将from元素替换为to元素
  x <- vector
  for (i in c(1:length(from))) {
    x[x == from[i]] <- to[i]
  }
  return(x)
}

inter <- function(list, ab, a, b, v) {      #计算交互效应，list列表，ab、a、b交互、因子1、因子2序号向量，v比较的数据
  out <- list()
  for (i in c(1:length(ab))) {
    up <- max(quantile(list[[a[i]]][v], 0.975, na.rm=T), quantile(list[[b[i]]][v], 0.975, na.rm=T))
    down <- min(quantile(list[[a[i]]][v], 0.025, na.rm=T), quantile(list[[b[i]]][v], 0.025, na.rm=T))
    E <- as.vector(unlist(list[[a[i]]][v] + list[[b[i]]][v]))
    Emean <- mean(E, na.rm=T)
    Edown <- quantile (E, 0.025, na.rm=T)
    Eup   <- quantile (E, 0.975, na.rm=T)
    list[[ab[i]]]$inter[unlist(list[[ab[i]]][v]) - E > Edown - Emean & unlist(list[[ab[i]]][v]) - E < Eup - Emean] <- "Additive"
    if (mean(unlist(list[[a[i]]][v]), na.rm=T) / mean(unlist(list[[b[i]]][v]), na.rm=T) < 0) {
      list[[ab[i]]]$inter[((unlist(list[[ab[i]]][v]) - E) - (Eup - Emean)) * Emean > 0 & 
                            (unlist(list[[ab[i]]][v]) - up) * (unlist(list[[ab[i]]][v]) - down) > 0] <- "+Synergistic"
      list[[ab[i]]]$inter[((unlist(list[[ab[i]]][v]) - E) - (Eup - Emean)) * Emean > 0 & 
                            (unlist(list[[ab[i]]][v]) - up) * (unlist(list[[ab[i]]][v]) - down) < 0] <- "-Antagonistic"
      list[[ab[i]]]$inter[((unlist(list[[ab[i]]][v]) - E) - (Eup - Emean)) * Emean < 0 & 
                            (unlist(list[[ab[i]]][v]) - up) * (unlist(list[[ab[i]]][v]) - down) < 0] <- "+Antagonistic"
      list[[ab[i]]]$inter[((unlist(list[[ab[i]]][v]) - E) - (Eup - Emean)) * Emean < 0 & 
                            (unlist(list[[ab[i]]][v]) - up) * (unlist(list[[ab[i]]][v]) - down) > 0] <- "-Synergistic"
    } else {
      list[[ab[i]]]$inter[list[[ab[i]]][v] - E < Edown - Emean] <- "Antagonistic"
      list[[ab[i]]]$inter[list[[ab[i]]][v] - E > Eup - Emean] <- "Synergistic"
    }
    out[[i]] <- list[[ab[i]]]
  }
  return(out)
}

si <- function(data, v, group) {               #将data数据框按group分组计算v的平均值mean，2.5%和97.5%的分位数down、up，和大于0的比例p
  library(plyr)                                #dlply拆分数据框
  datac <- ddply(data, group,
                 .fun = function(xx, col) {
                   c(mean = mean     (xx[[col]], na.rm=T),
                     down = quantile (xx[[col]], 0.025, na.rm=T),
                     up   = quantile (xx[[col]], 0.975, na.rm=T),
                     p    = length(which(xx[[col]]>0))/length(xx[[col]])
                   )
                 },
                 v
  )
  names(datac) <- c(group, "mean", "down", "up", "p")
  datac <- datac[which(is.nan(datac$mean) == F),]
  return(datac)
}

replace <- function(vector, from, to) {       #批量替换向量vector中的元素，将from元素替换为to元素
  x <- vector
  for (i in c(1:length(from))) {
    x[x == from[i]] <- to[i]
  }
  return(x)
}


##################
###  统计分析  ###
##################

#效应值计算：E-Effect size，P-Performance，R-RCI
##全球变化因子
E_P_ECF <- bootstrap(Data, "ReTG", "ReTC", group_c = "Target_Type", c("ECFs", "NUM_ECFs_Record"), 1000, 100)
E_C_ECF <- bootstrap(Data, "TGRCI", "TCRCI", group_c = "Target_Type", c("ECFs", "NUM_ECFs_Record"),  1000, 100)

#计算差异值Difference size：仅ECF-M与ECF-S之间
Index <- unique(rbindlist(E_P_ECF)[,c(1,3)])
Index$Num <- c(1:nrow(Index))
Index_M <- Index[nchar(Index$ECFs) == 2,]
for (i in c(1:nrow(Index_M))) {
  Index_M$ECF1[i] <- substring(Index_M[i,1], 1, 1)
  Index_M$ECF2[i] <- substring(Index_M[i,1], 2, 2)
}
Index_M$Num1[Index_M$Target_Type == "Invasive"] <- as.numeric(replace(Index_M$ECF1[Index_M$Target_Type == "Invasive"], c("C", "D", "N", "P", "T"), c(1, 3, 9, 13, 17)))
Index_M$Num2[Index_M$Target_Type == "Invasive"] <- as.numeric(replace(Index_M$ECF2[Index_M$Target_Type == "Invasive"], c("C", "D", "N", "P", "T"), c(1, 3, 9, 13, 17)))
Index_M$Num1[Index_M$Target_Type == "Native"] <- as.numeric(replace(Index_M$ECF1[Index_M$Target_Type == "Native"], c("C", "D", "N", "P", "T"), c(2, 4, 10, 14, 18)))
Index_M$Num2[Index_M$Target_Type == "Native"] <- as.numeric(replace(Index_M$ECF2[Index_M$Target_Type == "Native"], c("C", "D", "N", "P", "T"), c(2, 4, 10, 14, 18)))
D_P_ECF <- list()
for (i in c(1:nrow(Index_M))) {
  D_P_ECF[[i]] <- cbind(E_P_ECF[[Index_M$Num[i]]][,c(1:3)],
                        data.frame(D = E_P_ECF[[Index_M$Num[i]]]$E - 
                                     (E_P_ECF[[Index_M$Num1[i]]]$E + E_P_ECF[[Index_M$Num2[i]]]$E)))
}
D_C_ECF <- list()
for (i in c(1:nrow(Index_M))) {
  D_C_ECF[[i]] <- cbind(E_C_ECF[[Index_M$Num[i]]][,c(1:3)],
                        data.frame(D = E_C_ECF[[Index_M$Num[i]]]$E - 
                                     (E_C_ECF[[Index_M$Num1[i]]]$E + E_C_ECF[[Index_M$Num2[i]]]$E)))
}

#交互判定
I_P_ECF <- inter(E_P_ECF, 
                 c(5, 6, 7, 8, 11, 12, 15, 16, 19, 20, 21, 22, 23, 24), 
                 c(3, 4, 3, 4, 9, 10, 13, 14, 17, 18, 17, 18, 17, 18), 
                 c(1, 2, 9, 10, 1, 2, 9, 10, 1, 2, 3, 4, 9, 10), "E")
I_C_ECF <- inter(E_C_ECF, 
                 c(5, 6, 7, 8, 11, 12, 15, 16, 19, 20, 21, 22, 23, 24), 
                 c(3, 4, 3, 4, 9, 10, 13, 14, 17, 18, 17, 18, 17, 18), 
                 c(1, 2, 9, 10, 1, 2, 9, 10, 1, 2, 3, 4, 9, 10), "E")

#合并效应值、差异值
SE_P_ECF <- mq(rbindlist(E_P_ECF), "E", c("Target_Type", "ECFs", "NUM_ECFs_Record"))
SE_C_ECF <- mq(rbindlist(E_C_ECF), "E", c("Target_Type", "ECFs", "NUM_ECFs_Record"))
SD_P_ECF <- mq(rbindlist(D_P_ECF), "D", c("Target_Type", "ECFs"))
SD_C_ECF <- mq(rbindlist(D_C_ECF), "D", c("Target_Type", "ECFs"))
SI_P_ECF <- as.data.frame(table(rbindlist(I_P_ECF)$ECFs, rbindlist(I_P_ECF)$Target_Type, rbindlist(I_P_ECF)$inter))
names(SI_P_ECF) <- c("ECFs", "Target_Type", "inter", "Freq")
SI_C_ECF <- as.data.frame(table(rbindlist(I_C_ECF)$ECFs, rbindlist(I_C_ECF)$Target_Type, rbindlist(I_C_ECF)$inter))
names(SI_C_ECF) <- c("ECFs", "Target_Type", "inter", "Freq")
SE_P_F <- mq(rbindlist(E_P_ECF), "E", c("Target_Type", "NUM_ECFs_Record"))
SE_C_F <- mq(rbindlist(E_C_ECF), "E", c("Target_Type", "NUM_ECFs_Record"))
SD_P_F <- mq(rbindlist(D_P_ECF), "D", c("Target_Type", "NUM_ECFs_Record"))
SD_C_F <- mq(rbindlist(D_C_ECF), "D", c("Target_Type", "NUM_ECFs_Record"))

#获取样本量
N_P_ECF <- ddply(Data, c("Target_Type", "ECFs"),
                 .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_P_ECF)[3] <- "N"
N_C_ECF <- ddply(Data, c("Target_Type", "ECFs"),
                 .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_C_ECF)[3] <- "N"
N_P_F <- ddply(Data, c("Target_Type", "NUM_ECFs_Record"),
               .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_P_F)[3] <- "N"
N_C_F <- ddply(Data, c("Target_Type", "NUM_ECFs_Record"),
               .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_C_F)[3] <- "N"
N_P_N <- ddply(Data, c("Target_Type", "NFIX_Target"),
               .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_P_N)[3] <- "N"
N_C_N <- ddply(Data, c("Target_Type", "NFIX_Target"),
               .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_C_N)[3] <- "N"
N_P_LC <- ddply(Data, c("Target_Type", "Lifecycle_Target"),
                .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_P_LC)[3] <- "N"
N_C_LC <- ddply(Data, c("Target_Type", "Lifecycle_Target"),
                .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_C_LC)[3] <- "N"
N_P_FG <- ddply(Data, c("Target_Type", "FUN_Group_Target"),
                .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_P_FG)[3] <- "N"
N_C_FG <- ddply(Data, c("Target_Type", "FUN_Group_Target"),
                .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_C_FG)[3] <- "N"
N_P_I <- ddply(Data, c("Target_Type", "Type_IND"),
               .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_P_I)[3] <- "N"
N_C_I <- ddply(Data, c("Target_Type", "Type_IND"),
               .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_C_I)[3] <- "N"

#数据整理
P_F <- merge(N_P_F, SE_P_F, by = c("Target_Type", "NUM_ECFs_Record"), all = TRUE) 
C_F <- merge(N_C_F, SE_C_F, by = c("Target_Type", "NUM_ECFs_Record"), all = TRUE) 
P_F$ECFs[P_F$NUM_ECFs_Record == 1] <- "ECF(S)"
P_F$ECFs[P_F$NUM_ECFs_Record == 2] <- "ECF(M)"
C_F$ECFs[C_F$NUM_ECFs_Record == 1] <- "ECF(S)"
C_F$ECFs[C_F$NUM_ECFs_Record == 2] <- "ECF(M)"
P_ECF <- merge(N_P_ECF, SE_P_ECF, by = c("Target_Type", "ECFs"), all = TRUE) 
C_ECF <- merge(N_C_ECF, SE_C_ECF, by = c("Target_Type", "ECFs"), all = TRUE) 
P_ECF <- rbind(merge(N_P_ECF, SE_P_ECF, by = c("Target_Type", "ECFs"), all = TRUE), P_F)
C_ECF <- rbind(merge(N_C_ECF, SE_C_ECF, by = c("Target_Type", "ECFs"), all = TRUE), C_F)
names(SD_P_F)[2] <- "ECFs"
SD_P_F$ECFs[SD_P_F$ECFs == 2] <- "ECF(M)"
SD_P_ECF <- rbind(SD_P_ECF, SD_P_F)
names(SD_C_F)[2] <- "ECFs"
SD_C_F$ECFs[SD_C_F$ECFs == 2] <- "ECF(M)"
SD_C_ECF <- rbind(SD_C_ECF, SD_C_F)

#数据输出
##表格
write.xlsx(P_ECF, "1-result.xlsx", sheetName = "P_ECF")
write.xlsx(C_ECF, "1-result.xlsx", sheetName = "C_ECF", append = TRUE)
write.xlsx(SD_P_ECF, "1-result.xlsx", sheetName = "SD_P_ECF", append = TRUE)
write.xlsx(SD_C_ECF, "1-result.xlsx", sheetName = "SD_C_ECF", append = TRUE)
write.xlsx(SI_P_ECF, "1-result.xlsx", sheetName = "SI_P_ECF", append = TRUE)
write.xlsx(SI_C_ECF, "1-result.xlsx", sheetName = "SI_C_ECF", append = TRUE)



##绘图
###绘图设置
cc <- c("#00afbb", "#1bbf93", "#1bbf93", "#86c44e", "#e7b800", "#e7b800", "#e7b800", "#e7b800")
ce <- c("#5db5bb", "#66c5a3", "#66c5a3", "#9ecd7a", "#e7ca5d", "#e7ca5d", "#e7ca5d", "#e7ca5d")
ct <- c("#000000", "#000000", "#000000")
c3 <- c("#1f78b4", "#33a02c")
c7 <- c("#007f87", "#fF1.4f08", "#00afbb", "#e7b800", "#fc4e08", "#c93e06", "#00d2e1")
#c7 <- c("#1f78b4", "#fF1.4f08", "#00afbb", "#e7b800", "#fc4e08", "#c93e06", "#6a3d9a")
windowsFonts(myFont = windowsFont("Times New Roman")) #字体修改

###制图
####图1-3单因子结果
D1.3 <- rbind(
  cbind(P_ECF[P_ECF$NUM_ECFs_Record == 1,], 
        data.frame(PC = rep("(A) Performance", length(P_ECF[P_ECF$NUM_ECFs_Record == 1,1])))),
  cbind(C_ECF[C_ECF$NUM_ECFs_Record == 1,], 
        data.frame(PC = rep("(B) Competitiveness", length(C_ECF[C_ECF$NUM_ECFs_Record == 1,1]))))
)
D1.3$PC <- factor(D1.3$PC, levels = c("(A) Performance", "(B) Competitiveness"))
D1.3$sig[D1.3$p < 0.025 | D1.3$p > 0.975] <- "*"

F1.3 <- ggplot(data = D1.3)+
  geom_point(aes(x = mean, y = ECFs, color = Target_Type, shape = Target_Type),
             size = 3, position = position_dodge(width = -0.5))+
  geom_errorbarh(aes(y = ECFs, color = Target_Type, xmin = down, xmax = up), 
                 height = 0.2, size = 0.8, position = position_dodge(width = -0.5))+
  geom_vline(xintercept = 0, linetype = "dashed", size = 0.7)+
  geom_text(aes(x = up, y = ECFs, color = Target_Type, label = N),
            position=position_dodge(width = -0.5), hjust = -1, show.legend = FALSE)+
  geom_text(aes(x = down, y = ECFs, color = Target_Type, label = sig), size = 8,
            position=position_dodge(width = -0.5), hjust = 1.5, vjust = 0.7, show.legend = FALSE)+
  scale_x_continuous(expand = c(0.1,0.05))+
  scale_color_manual(values = c3)+
  scale_fill_manual(values = c3)+
  theme_bw()+
  scale_y_discrete(limits = factor(c("ECF(S)","C","N","D","P","T"), 
                                   levels=c("ECF(S)","C","N","D","P","T")))+
  labs(x = "Effect size")+
  facet_wrap(vars(PC),scales= "free_x")+
  theme(axis.title.y=element_blank(),
        legend.background=element_blank(),
        legend.title=element_blank(),
        legend.text=element_text(family = "myFont", size = 12),
        legend.justification=c(0,0), 
        legend.position=c(0.01,0.1),
        title=element_text(family = "myFont",size=14,),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12),
        strip.text.x = element_text(family = "myFont",size=14))

####图1.4两因子结果
D1.4 <- rbind(
  cbind(P_ECF[P_ECF$NUM_ECFs_Record == 2,], 
        data.frame(PC = rep("(A) Performance", length(P_ECF[P_ECF$NUM_ECFs_Record == 2,1])))),
  cbind(C_ECF[C_ECF$NUM_ECFs_Record == 2,], 
        data.frame(PC = rep("(B) Competitiveness", length(C_ECF[C_ECF$NUM_ECFs_Record == 2,1]))))
)
D1.4$PC <- factor(D1.4$PC, levels = c("(A) Performance", "(B) Competitiveness"))
D1.4$sig[D1.4$p < 0.025 | D1.4$p > 0.975] <- "*"

F1.4 <- ggplot(data = D1.4)+
  geom_point(aes(x = mean, y = ECFs, color = Target_Type, shape = Target_Type),
             size = 3, position = position_dodge(width = -0.5))+
  geom_errorbarh(aes(y = ECFs, color = Target_Type, xmin = down, xmax = up), 
                 height = 0.2, size = 0.8, position = position_dodge(width = -0.5))+
  geom_vline(xintercept = 0, linetype = "dashed", size = 0.7)+
  geom_text(aes(x = up, y = ECFs, color = Target_Type, label = N),
            position=position_dodge(width = -0.5), hjust = -1, show.legend = FALSE)+
  geom_text(aes(x = down, y = ECFs, color = Target_Type, label = sig), size = 8,
            position=position_dodge(width = -0.5), hjust = 1.5, vjust = 0.7, show.legend = FALSE)+
  scale_x_continuous(expand = c(0.1,0.05))+
  scale_color_manual(values = c3)+
  scale_fill_manual(values = c3)+
  theme_bw()+
  scale_y_discrete(limits = factor(c("ECF(M)","NC","DC","DN","PN","TC","TN","TD"), 
                                   levels=c("ECF(M)","NC","DC","DN","PN","TC","TN","TD")))+
  labs(x = "Effect size")+
  facet_wrap(vars(PC),scales= "free_x")+
  theme(axis.title.y=element_blank(),
        legend.background=element_blank(),
        legend.title=element_blank(),
        legend.text=element_text(family = "myFont", size = 12),
        legend.justification=c(0,0), 
        legend.position=c(0.01,0.1),
        title=element_text(family = "myFont",size=14,),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12),
        strip.text.x = element_text(family = "myFont",size=14))

####图1.5交互结果
D1.5 <- rbind(
  cbind(SD_P_ECF, 
        data.frame(PC = rep("(A) Performance", length(P_ECF[P_ECF$NUM_ECFs_Record == 2,1])))),
  cbind(SD_C_ECF, 
        data.frame(PC = rep("(B) Competitiveness", length(C_ECF[C_ECF$NUM_ECFs_Record == 2,1]))))
)
D1.5$PC <- factor(D1.5$PC, levels = c("(A) Performance", "(B) Competitiveness"))
D1.5$sig[D1.5$p < 0.025 | D1.5$p > 0.975] <- "*"
D1.5$N <- D1.4$N

F1.5 <- ggplot(data = D1.5)+
  geom_point(aes(x = mean, y = ECFs, color = Target_Type, shape = Target_Type),
             size = 3, position = position_dodge(width = -0.5))+
  geom_errorbarh(aes(y = ECFs, color = Target_Type, xmin = down, xmax = up), 
                 height = 0.2, size = 0.8, position = position_dodge(width = -0.5))+
  geom_vline(xintercept = 0, linetype = "dashed", size = 0.7)+
  geom_text(aes(x = up, y = ECFs, color = Target_Type, label = N),
            position=position_dodge(width = -0.5), hjust = -1, show.legend = FALSE)+
  geom_text(aes(x = down, y = ECFs, color = Target_Type, label = sig), size = 8,
            position=position_dodge(width = -0.5), hjust = 1.5, vjust = 0.7, show.legend = FALSE)+
  scale_x_continuous(expand = c(0.1,0.05))+
  scale_color_manual(values = c3)+
  scale_fill_manual(values = c3)+
  theme_bw()+
  scale_y_discrete(limits = factor(c("ECF(M)","NC","DC","DN","PN","TC","TN","TD"), 
                                   levels=c("ECF(M)","NC","DC","DN","PN","TC","TN","TD")))+
  labs(x = "Different size")+
  facet_wrap(vars(PC),scales= "free_x")+
  theme(axis.title.y=element_blank(),
        legend.background=element_blank(),
        legend.title=element_blank(),
        legend.text=element_text(family = "myFont", size = 12),
        legend.justification=c(0,0), 
        legend.position=c(0.01,0.1),
        title=element_text(family = "myFont",size=14,),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12),
        strip.text.x = element_text(family = "myFont",size=14))



###存图
ggsave(file = "F1.3.svg",plot = F1.3, width = 8, height = 6)
ggsave(file = "F1.4.svg",plot = F1.4, width = 8, height = 7.5)
ggsave(file = "F1.5.svg",plot = F1.5, width = 8, height = 7.5)

ggsave(file = "F1.3.tiff",plot = F1.3, width = 8, height = 6, bg = "white")
ggsave(file = "F1.4.tiff",plot = F1.4, width = 8, height = 7.5, bg = "white")
ggsave(file = "F1.5.tiff",plot = F1.5, width = 8, height = 7.5, bg = "white")

